
MODULE MOD_SOLAR
# if defined (HEATING_SOLAR)
  USE ALL_VARS
  USE MOD_UTILS
  USE MOD_PREC

  IMPLICIT NONE
  SAVE

  LOGICAL HEATING_SOLAR_ON
  CHARACTER(LEN=80) HEATING_SOLAR_TYPE
  CHARACTER(LEN=80) HEATING_SOLAR_FILE
  CHARACTER(LEN=80) HEATING_SOLAR_KIND
  REAL(SP) :: ZM                            !!ZM - ANEMOMETER HEIGHT (M)
  REAL(SP) :: LATITUDE_REFERENCE
  REAL(SP) :: LONGITUDE_REFERENCE
  INTEGER  :: JULIANDAY_REFERENCE   
  NAMELIST /NML_HEATING_SOLAR/        &
       & HEATING_SOLAR_ON,            &
       & HEATING_SOLAR_TYPE,          &
       & HEATING_SOLAR_FILE,          &
       & HEATING_SOLAR_KIND,          &
       & ZM,                          &
       & LATITUDE_REFERENCE,          &
       & LONGITUDE_REFERENCE,         &
       & JULIANDAY_REFERENCE   
   
   REAL(SP), ALLOCATABLE :: CORRG(:),CORR(:)           !!LATITUDE OF NODES
    
  CONTAINS

!==============================================================================|
!   Input Parameters Which Control the Calculation of Heat Flux                |
!==============================================================================|

   SUBROUTINE HEATING_SOLAR_NAMELIST_INITIALIZE
   USE control, only : casename
   IMPLICIT NONE
   
   HEATING_SOLAR_ON    = .FALSE.                
   HEATING_SOLAR_TYPE  = "'flux' or 'body'"              
   HEATING_SOLAR_FILE  = trim(casename)//"_hfx.nc"              
   HEATING_SOLAR_KIND  = "Options:"//TRIM(CNSTNT)//","//TRIM(STTC)//","//TRIM(TMDPNDNT)//","//TRIM(PRDC)//","//TRIM(VRBL)              
   ZM                  = 10.  ! Unit = m                                 
   LATITUDE_REFERENCE  = 0.0_SP
   LONGITUDE_REFERENCE = 0.0_SP
   JULIANDAY_REFERENCE = 0   
   
   RETURN
   END SUBROUTINE HEATING_SOLAR_NAMELIST_INITIALIZE
   
!------------------------------------------------------------------------------|
   SUBROUTINE HEATING_SOLAR_NAMELIST_PRINT
   USE CONTROL, ONLY : IPT
   
   IMPLICIT NONE
   
   WRITE(UNIT=IPT,NML=NML_HEATING_SOLAR)
   
   RETURN
   END SUBROUTINE HEATING_SOLAR_NAMELIST_PRINT  
   
!------------------------------------------------------------------------------|
   SUBROUTINE HEATING_SOLAR_NAMELIST_READ    
   USE CONTROL, ONLY : casename,NMLUNIT
   
   IMPLICIT NONE
   
   INTEGER :: IOS, I
   CHARACTER(LEN=120) :: FNAME
   CHARACTER(LEN=160) :: PATHNFILE
   
   IF(DBG_SET(DBG_SBR)) &
         & WRITE(IPT,*) "Subroutine Begins: Read_Heating_Solar_Namelist;"

   IOS = 0

   FNAME = "./"//trim(casename)//"_run.nml"

   CALL FOPEN(NMLUNIT,trim(FNAME),'cfr')

   !READ NAME LIST FILE
    REWIND(NMLUNIT)

   ! Read IO Information
   READ(UNIT=NMLUNIT, NML=NML_HEATING_SOLAR,IOSTAT=ios)
   if(ios .NE. 0 ) Then
     if(DBG_SET(dbg_log)) write(UNIT=IPT,NML=NML_HEATING_SOLAR)
     Call Fatal_error("Can Not Read NameList NML_HEATING_SOLAR from file: "//trim(FNAME))
   end if
   CLOSE(NMLUNIT)
   
   IF(HEATING_SOLAR_ON .AND. .NOT. HEATING_ON)THEN
     HEATING_ON = .TRUE.
   ELSE IF(HEATING_SOLAR_ON .AND. HEATING_ON)THEN
     CALL FATAL_ERROR("IN NAMELIST CONTROL FILE, IF HEATING_SOLARCALCULATE_ON = TRUE, ", &
                      "THEN SET HEATING_ON = FALSE") 
   END IF 
   IF(HEATING_SOLAR_TYPE == 'flux') HEATING_TYPE = 'flux'
   IF(HEATING_SOLAR_TYPE == 'body') HEATING_TYPE = 'body'

   RETURN
   END SUBROUTINE HEATING_SOLAR_NAMELIST_READ

!==============================================================================|

!---------------------------------------------------------------------------
      SUBROUTINE SOLAR(TA,DP,WS,CC,TWAT,CH,CD,RLATD,RLONG,IJD0,MONTH,HR,SURFLX,Q0,HS,HL,RLN)
            

!
!++++++++++++++++++++++
!
!     SUBROUTINE SOLAR
!
!++++++++++++++++++++++
!
!
!. PURPOSE:    TO CALCULATE THE INSOLATION AND SURFACE HEAT FLUX
!              GIVEN THE AIR TEMPERATURE, WIND SPEED, DEW POINT OR
!              VAPOR PRESSURE, WATER TEMPERATURE, DRAG COEFFICIENTS, 
!              JULIAN DATE, AND CLOUD COVERAGE
!
!. ALGORITHMS: SENSIBLE AND LATENT HEAT FLUXES ARE CALCULATED FROM
!              BULK TRANSFER EQUATIONS.  LONGWAVE RADIATION LOSS
!              IS FROM WYRTKI (1965).
!              GLOBAL RADIATION IS FROM GUTTMAN AND MATTHEWS
!              (1979), IVANOFF (1977) AND COTTON (1979). 
!              GLOBAL RAD IS FROM 290 - 4000  
!              MMICRONS.  SOLAR ALBEDO IS FROM AUBERT AND RICHARDS (1981)
!              (P. 173).
!       
!. HISTORY:    ORIGINAL WRITTEN BY M.J.McCORMICK OF GLERL
!              ; MODIFIED FOR USE WITH GLFS BY D.WELSH OF OHIO STATE,1990
!              ; MODIFIED BY D.WELSH 4/93 FOR USE WITH NEW VERSION OF
!                  CIRCULATION MODEL. JD & IHR NOW IN CALL (CALLED BY
!                  TAU) ; MONTH NOW CALCULATED IN SOLAR ; LONGITUDE &
!                  LATITUDE NOW SET = PRINCETON MODEL'S ALON & ALAT ARRAYS.
!                 - THESE CHANGES / ADDITIONS ARE INDICATED THUS :
!
!                   C++++++++++++++++++++++++++++++++
!
!                        CHANGES/ADDITIONS HERE
!
!         4/97 DJS:   ADDED Q0,RLN,HS,AND HL AS RETURNED PARAMETERS
!         10/15/2007 DJS: NEW VERSION FOR USE IN FVCOM
!                   C++++++++++++++++++++++++++++++++
!  
!. PARAMETERS:(ALL UNITS ARE MKS)
!
!   CALLING VARIABLES (*=INPUT):
!    *TA - AIR TEMPERATURE
!    *DP - DEW POINT
!    *WS - WIND SPEED
!    *CC - CLOUD COVER (0-1)
!    *TWAT - WATER TEMPERATURE
!    *CH - BULK AERODYNAMIC COEFFICIENT FOR HEAT
!    *CD - BULK AERODYNAMIC COEFFICIENT FOR MOMENTUM
!    *RLATD - LATITUDE
!    *RLONG - LONGITUDE (WEST OF GREENWICH)
!    *IJD0 - DAY OF YEAR
!    *HR - HOUR OF DAY (LOCAL STANDARD TIME)
!    SURFLXC - TOTAL HEAT FLOUX (W/M**2)
!    Q0 - SHORT WAVE RADIATION (W/M**2)
!
!..INTERNAL VARIABLES:
!.. CP - SPECIFIC HEAT OF DRY AIR AT CONSTANT PRESSURE (J/(KG C))
!... CW - SPECIFIC HEAT OF WATER (J/(KG C))
!..  CP AND CW FROM BRUTSAERT (1982)
!... SIGMA - STEFAN-BOLZMANN CONSTANT (W/M**2)
!... SOLC - SOLAR CONSTANT (W/M**2) (FROHLICH, 1977)
!... TRANS - AVERAGE ANNUAL SOLAR TRANSMISSION (DECIMAL), (REMNANT
!...        THINKING, NO LONGER USED)
!... EXT1 - EXTINCTION COEFFICIENT FOR VISIBLE LIGHT IN WATER (1/M)
!... EXT2 - EXTINCTION COEFFICIENT FOR IR RADIATION IN WATER (1/M)
!..  SURFLX - SURFACE HEAT FLUX
!..  WS - WIND SPEED (M/S)
!..  TIME - HOURS SINCE BEGINNING OF FIRST YEAR
!..  DP - DEW POINT (DEG C)
!..  SC - CLOUD AMOUNT COVERING SKY (DECIMAL)
!..  RHOA - AIR DENSITY (KG/M**3)
!..  RHOW(T,Z) - WATER DENSITY AS A FUNCTION OF T, AND Z (KG/M**3)
!..  CH - BULK AERODYNAMIC COEFFICIENT FOR HEAT
!..  CD - "       "          "          "  MOMENTUM
!..  RAD - MEASURED GLOBAL RADIATION (W/M**2)
!..  MISS - NUMBER OF TIME STEPS RAD DATA ARE MISSING
!..  DEW - LOGICAL VARIABLE. IF .TRUE. THEN HUMIDITY UNITS ARE
!..        MBAR, OTHERWISE, DEW POINT TEMPERATURE.
!..  ES - ELLIPTICITY OF EARTH'S ORBIT
!..  ED - 1-(SQUARE OF ES)
!..  RM - TRUE SOLAR NOON
!..  SRCS - CLEAR SKY SOLAR RADIATION (COTTON, 1979)
!..  CLD - PERCENT TRANSMISSION OF SRCS DUE TO CLOUDS
!..  C'S ANS A'S - COEFFICIENTS DERIVED BY COTTON (1979)
!..                FOR USE IN HIS REGRESSION RELATION.
!.. N.B.  THESE JUST MENTIONED COEFFICIENTS ARE FOR MADISON,
!.. WISC. AND ARE USEFUL FOR HOURLY RADIATION (GLOBAL) ESTIMATES.
!..................................................................

  IMPLICIT NONE
      REAL(SP), INTENT(IN) :: TA,DP,WS,CC,TWAT,CH,CD,RLATD,RLONG,HR
      INTEGER, INTENT(IN) :: IJD0,MONTH
      REAL(SP), INTENT(OUT) :: SURFLX,Q0,HL,HS,RLN
      LOGICAL DEW
      REAL(SP) :: CP, CW, SIGMA
      REAL(SP) :: SOLC,PHASE,TWOPI
      REAL(SP) :: DEG
      REAL(SP) :: ES, ED
      INTEGER,DIMENSION(12) :: IMON
      REAL(SP),DIMENSION(12,2) :: A0
      REAL(SP),DIMENSION(2) :: A1,A2,A3
      REAL(SP) :: C0,C1,C2,C3
      REAL(SP) :: SC,SST,TAK,DPK,TWK,RHLAT,RHOA,ETW,EDP,QW,QA,GMT,DAY
      REAL(SP) :: ANGD,SIG,STHETA,ARCS,CTHETA,RM,HRANG,SINLAT,COSLAT,CST,SOL
      REAL(SP) :: SRCS,CLD,ALB,ALBEDO
      INTEGER :: JD,MONTHR,M,N
!
      DATA CP, CW, SIGMA/1005., 4186., 5.6697E-08/
      DATA SOLC,PHASE,TWOPI/1375., .4583, 6.28319/
      DATA DEG/273.15/
      DATA ES, ED/0.0167238, .99986/
      DATA IMON/32,60,91,121,152,182,213,244,274,305,335,366/
      DATA A0/10.,37.,62.,-6.,-84.,-157.,-180.,-164.,-103.,-61.,-12.,6.,-18.,10.,12.,-56.,-148.,-217.,-222.,-216.,-178.,-106.,-59.,-50./
      DATA A1/2460., 2550./
      DATA A2/3010., 2800./
      DATA A3/-1750., -1580./
      DATA C0,C1,C2,C3/.999, -.425, .922, -1.14/
      SC=CC
!
      DEW=.FALSE.
      TWOPI=6.28318531
      SST=TWAT
      TAK=TA + DEG
      DPK=DP + DEG
      TWK=SST + DEG
!      JDR=IJD0
!      SHR=HR
!      WRITE(6,*)'IJD0 : ',IJD0
!      WRITE(6,*)'HR : ',HR
!.
!.. CALCULATE LATENT HEAT OF VAPORIZATION
!.
      RHLAT=(2.501-.00236*SST)*1.E+06
      RHOA=354.65/TAK
      ETW=VAP(TWK)
      IF(DEW) THEN
      EDP=DP
      ELSE
      EDP=VAP(DPK)
      END IF 
!.
!.. CALCULATE SPECIFIC HUMIDITY AT WATER SURFACE AND AT INSTRUMENT
!.. HEIGHT
!.
      QW=.622*ETW/(1018.-.378*ETW)
      QA=.622*EDP/(1018.-.378*EDP)
!.
!.. CALCULATE LATENT HEAT LOSS
!.
      HL=-CD*RHOA*RHLAT*WS*(QW-QA)
!.
!.. CALCULATE SENSIBLE HEAT LOSS
!.
      HS=-CH*CP*RHOA*WS*(SST-TA)
!.
!.. CALCULATE LONGWAVE RADIATION LOSS
!.
      RLN=-SIGMA*TWK*TWK*TWK*TWK*(0.39-.05*SQRT(VAP(TAK)))*(1.-.67*SC*SC) - 4.*SIGMA*TWK*TWK*TWK*(SST-TA)
!.
!.. SHORT WAVE GLOBAL RADIATION
!.. CALCULATE COSINE OF ZENITH ANGLE, CST
!.. THETA = SUN'S DECLINATION, HR = LOCAL HOUR ANGLE
!.
      JD=IJD0
!      JDRR=JD
!      DHOUR=HR
      GMT=HR
!      DHRR=GMT
!.
      DAY=FLOAT(JD)
!      DAYR=DAY
!JQI <
!      IF(JD.LE.31)MONTH=1
!      IF((JD.GT.31).AND.(JD.LE.59))MONTH=2
!      IF((JD.GT.59).AND.(JD.LE.90))MONTH=3
!      IF((JD.GT.90).AND.(JD.LE.120))MONTH=4
!      IF((JD.GT.120).AND.(JD.LE.151))MONTH=5
!      IF((JD.GT.151).AND.(JD.LE.181))MONTH=6
!      IF((JD.GT.181).AND.(JD.LE.212))MONTH=7
!      IF((JD.GT.212).AND.(JD.LE.243))MONTH=8
!      IF((JD.GT.243).AND.(JD.LE.273))MONTH=9
!      IF((JD.GT.273).AND.(JD.LE.304))MONTH=10
!      IF((JD.GT.304).AND.(JD.LE.334))MONTH=11
!      IF(JD.GT.334)MONTH=12
!JQI >
      MONTHR=MONTH
!
!..CHANGED FROM VAX VERSION SO THAT ANGLES ARE NOW IN RADIANS
!
      ANGD=TWOPI*(DAY-1.)/365.242
!      ANGDR=ANGD
      SIG=279.9348+1.914827*SIN(ANGD) - 0.079525*COS(ANGD)+0.01993*SIN(2.*ANGD)-0.00162*COS(2.*ANGD) + ANGD*360./TWOPI
!      SIGR=SIG
      SIG=SIG*TWOPI/360.
      STHETA=SIN(0.4091722)*SIN(SIG)
!      STHETAR=STHETA
      ARCS=ASIN(STHETA)
!      ARCSR=ARCS
      CTHETA=COS(ARCS)
!      CTHETAR=CTHETA
      RM=12.+(0.12357*SIN(ANGD)-0.004289*COS(ANGD)+0.153809*SIN(2.*ANGD)+.060783*COS(2.*ANGD))
!      RMR=RM
      HRANG=15.*(GMT-RM) - RLONG
!      HRANGR=HRANG
      HRANG=HRANG*TWOPI/360.
      SINLAT=SIN(RLATD*TWOPI/360.0)
!      SINLATR=SINLAT
      COSLAT=COS(RLATD*TWOPI/360.0)
!      COSLATR=COSLAT
      CST=SINLAT*STHETA+COSLAT*CTHETA*COS(HRANG)
!      CSTR=CST
      SOL=SOLC*((1.+ES*COS(ANGD))/ED)**2
!      SOLR=SOL
!
      Q0=0.
      SRCS=0.0
      CLD=0.0
!.
!.. IF SUN IS BELOW HORIZON OR ALBEDO IS GREATER THAN 1, Q = 0.
!.
!.. CHECK IF MEASURED OR MODELED RADIATION DATA ARE AVAILABLE.
!.. IF DATA ARE AVAILABLE USE IN PLACE OF PRESENT ALGORITHM.
!.
      IF(CST .GT. 0.) THEN
        ALB=.045/CST
        ALBEDO=ALB
      !EJA - fixed assumed bug in ALBEDO - 05/03/2016
      !ALBEDO=0.
        IF(ALB .GT. 1.) ALBEDO=1.
!.  
        CLD=C0 +C1*SC +C2*SC*SC +C3*SC*SC*SC
        M=MONTH
        N=1
        IF(RM .GT. 12.) N=2
        SRCS=A0(M,N) +A1(N)*CST +A2(N)*CST*CST +A3(N)*CST*CST*CST
!
        IF(SRCS.LT.0.0)SRCS=0.0
!.
!..CONVERT FROM KJ/HR TO W/M**2 AND INCLUDE ALBEDO
!.
        Q0=SRCS*CLD*(1.-ALBEDO)/3.6
!      MR=MONTH
!      NR=N
!      CLDR=CLD
!      ALBEDOR=ALBEDO
!      SRCSR=SRCS
!.
      ENDIF
!     IF(RAD .NE. -999.) Q0=RAD
!.
!.
      SURFLX=HS+HL+RLN+Q0

      RETURN
      END SUBROUTINE SOLAR

      REAL FUNCTION VAP(T)
!.
!.. FUNCTION CALCULATES SATURATION VAPOR PRESSURE IN MB.
!.. ALGORITHM FROM FLEAGLE AND BUSINGER (1980) P. 72.
!.  'T' IS IN DEGREES KELVIN
!.
     IMPLICIT NONE
     REAL(SP) :: T

     VAP=.01*EXP(26.25 - 5418./T)
!.
     RETURN
     END FUNCTION VAP

!--------------------------------------------------------------------
     SUBROUTINE UZL10(WS, TD, CD, CH)
!
! PURPOSE:
!               TO CALCULATE THE BULK AERODYNAMIC COEFFICIENTS FOR
!               MOMENTUM AND HEAT OVER A LAKE SURFACE AS FUNCTIONS
!               OF 10M WIND SPEED AND AIR-SEA TEMPERATURE DIFFERENCE.
!
! ALGORITHM:
!               TABLES OF CD AND CH FOR WS FROM 0
!               TO 50 M/S AND TD FROM -50 TO 50 C ARE CALCULATED
!               THE FIRST TIME UZL10 IS CALLED.  IN SUBSEQUENT
!               CALLS, CD AND CH ARE CALCULATED FROM THE TABLE
!               USING BILINEAR INTERPOLATION.
!
! ARGUMENTS:
!  ON INPUT:
!               WS - WIND SPEED (M / S)
!               TD - AIR-SEA TEMPERATURE DIFFERENCE (DEG K)
!  ON OUTPUT:
!               CD - BULK AERODYNAMIC COEFFICIENT FOR MOMENTUM
!               CH - BULK AERODYNAMIC COEFFICIENT FOR HEAT
!
! HISTORY:
!               WRITTEN BY D.J.SCHWAB, 1993, GLERL, ANN, ARBOR, MI
      IMPLICIT NONE
      SAVE

      REAL(SP),INTENT(IN) :: WS, TD
      REAL(SP),INTENT(OUT) :: CD, CH

      REAL(SP) CD10(51,101),CH10(51,101)
      INTEGER :: INIT,I,J
      REAL(SP) :: Z,WSX,TDX,Z0,FL,WI,WJ

      DATA INIT /0/
      IF(INIT.EQ.0) THEN
        INIT=1
        Z=10
        DO I=1,51
          WSX=FLOAT(I-1)
          DO J=1,101
            TDX=FLOAT(J-51)
            CALL UZL(WSX, Z, TDX, CD, CH, Z0, FL)
            CD10(I,J)=CD
            CH10(I,J)=CH
          ENDDO
        ENDDO
      ENDIF
      IF(WS.LT.0.) THEN
       PRINT *,'***WS LT 0 IN SUBROUTINE UZL10***'
       CD=0
       CH=0
       RETURN
      ENDIF
!cwpo      IF(ABS(TD.GT.50.)) THEN
      IF(ABS(TD).GT.50.) THEN
       PRINT *,'***ABS(TD) GT 50 IN SUBROUTINE UZL10***'
       CD=0
       CH=0
       RETURN
      ENDIF
      I=INT(WS)+1
      WI=WS-I+1
      J=INT(TD+50)+1
      WJ=(TD+50)-J+1
      CD=(1.-WI)*(1.-WJ)*CD10(I,J)+WI*(1.-WJ)*CD10(I+1,J)+(1.-WI)*WJ*CD10(I,J+1)+WI*WJ*CD10(I+1,J+1)
      CH=(1.-WI)*(1.-WJ)*CH10(I,J)+WI*(1.-WJ)*CH10(I+1,J)+(1.-WI)*WJ*CH10(I,J+1)+WI*WJ*CH10(I+1,J+1)
      IF(CD.GT.1.) PRINT *,'XXX',I,J,WI,WJ,CD,CD10(I,J),CD10(I+1,J),CD10(I,J+1),CD10(I+1,J+1),WS,TD
      RETURN
      END SUBROUTINE UZL10

      SUBROUTINE UZL(U, ZM, TD, CD, CH, Z0, FL)
!
! PURPOSE:
!               TO CALCULATE THE BULK AERODYNAMIC COEFFICIENTS FOR
!               MOMENTUM AND HEAT OVER A LAKE SURFACE AS FUNCTIONS
!               OF WIND SPEED AND AIR-SEA TEMPERATURE DIFFERENCE.
!
! ALGORITHM:
!               THERE IS AN OUTER ITERATION IN WHICH THE ROUGHNESS
!               LENGTH IS VARIED ACCORDING TO CHARNOCK'S FORMULA AND
!               AN INNER ITERATION IN WHICH THE STABILITY LENGTH
!               (MONIN-OBUKHOV LENGTH) IS VARIED ACCORDING TO THE
!               BUSINGER-DYER FORMULATION.  THE CONSTANT IN CHARMOCK'S
!               FORMULA IS CHOSEN SO THAT UNDER NEUTRAL CONDITIONS THE
!               10 M DRAG COEFFICIENT IS 0.0016.
!
! ARGUMENTS:
!  ON INPUT:
!               U - WIND SPEED (M / S)
!               ZM - ANEMOMETER HEIGHT (M)
!               TD - AIR-SEA TEMPERATURE DIFFERENCE (DEG K)
!  ON OUTPUT:
!               CD - BULK AERODYNAMIC COEFFICIENT FOR MOMENTUM
!               CH - BULK AERODYNAMIC COEFFICIENT FOR HEAT
!               Z0 - ROUGHNESS LENGTH (M)
!               FL - STABILITY LENGTH (M)
!
! HISTORY:
!               WRITTEN BY J. R. BENNETT AND J. D. BOYD, 1979, GLERL,
!               ANN, ARBOR, MI;  BASED ON A CONSTANT ROUGHNESS VERSION
!               WRITTEN BY PAUL LONG AND WILL SHAFFER OF THE TECHNIQUES
!               DEVELOPMENT LABORATORY, NOAA, SILVER SPRINGS, MD.
!
!               MODIFIED WHEN TRANSFERRED TO VAX FROM CDC IN 1984 TO WORK
!                CORRECTLY WITH 32 BIT REALS
!
!               MODIFIED 1/86 TO PREVENT MODIFICATION OF FIRST ARGUMENT IN
!                CALLING PARAMETERS (U) WHEN U < 0.001

      IMPLICIT NONE
      REAL(SP),INTENT(IN) :: U, ZM, TD
      REAL(SP),INTENT(OUT) :: CD, CH, Z0, FL
      REAL(SP) :: C1, C2, C3
      REAL(SP) :: B1, B2, B3
      REAL(SP) :: EPS,UM,FK,TBAR,ALPHA,BETA,GAMM,GAMT,UST1,H,DTHETA
      REAL(SP) :: S,X,FLI,X1,ARG1,ARG2,A,B,AA,BB,CC,ROOT,X2,TSTAR,USTAR,Z0NEW
      INTEGER :: ITER,I
      LOGICAL :: flag_30

      DATA C1, C2, C3 /.684E-4, 4.28E-3, -4.43E-4/
      DATA B1, B2, B3 /1.7989E-3, 4.865E-4, 3.9028E-5/
      EPS = .01
      UM = U
      IF (UM .LT. .001) UM = .001
      FK = .35
      TBAR = 278.
      ALPHA = 4.7
      BETA = .74
      GAMM = 15.
      GAMT = 9.
      UST1 = 0.04 * UM
      H = ZM
      DTHETA = TD
      IF (ABS(DTHETA) .LT. 1.E-7) DTHETA = SIGN(1.E-7,DTHETA)
!
! INITIAL GUESS FOR Z0
!
      Z0 = .00459 * UST1 * UST1
      S = UM * UM * TBAR / (9.8*DTHETA)
      IF (ABS(S) .GT. 1.E6) S = SIGN(1.E6,S)
# if defined (DOUBLE_PRECISION)
      X = DLOG(H/Z0)
# else
      X = ALOG(H/Z0)
# endif      
!
!       INITIAL GUESS FOR L
!
      FL = S / X
      DO ITER = 1, 20
        flag_30=.FALSE.
# if defined (DOUBLE_PRECISION)
        X = DLOG(H/Z0)
# else
        X = ALOG(H/Z0)
# endif      
        IF (ABS(FL) .GT. 500.) FL = SIGN(500.,FL)
        IF (FL .LE. 0.) THEN
!
! UNSTABLE SECTION (L LT 0 OR DT LT 0)
!
        FLI = 1. / FL
!
! ASSUME 5 ITERATIONS SUFFICIENT
!
          DO I = 1, 5
            X1 = GAMT * FLI
            ARG1 = SQRT(1. - X1*H)
            ARG2 = SQRT(1. - X1*Z0)
            A = BETA * LOG((ARG1 - 1.)*(ARG2 + 1.)/((ARG1 + 1.)*(ARG2 -0.999999)))
            X1 = GAMM * FLI
            ARG1 = (1. - X1*H) ** (.25)
            ARG2 = (1. - X1*Z0) ** (.25)
            B = LOG((ARG1 - 1.)*(ARG2 + 1.)/((ARG1 + 1.)*(ARG2 - 0.999999))) +2. * (ATAN(ARG1) - ATAN(ARG2))
            FL = S * A / (B*B)
            IF (ABS(FL) .GT. 500.) FL = SIGN(500.,FL)
            FLI = 1. / FL
          ENDDO
        ELSE
!
! STABLE SECTION
!
! TRY MILDLY STABLE-
!
          DO WHILE(.TRUE.)
            AA = X * X
            X1 = H - Z0
            BB = 9.4 * X1 * X - .74 * S * X
            CC = 4.7 * X1
            CC = CC * CC - CC * S
            ROOT = BB * BB - 4. * AA * CC
            IF (ROOT .LT. 0.) THEN
              flag_30=.TRUE.
              EXIT
            ENDIF
            FL = (-BB + SQRT(ROOT)) / (2.*AA)
            IF (FL .LE. H) THEN
              flag_30=.TRUE.
              EXIT
            ENDIF
            B = X + 4.7 * X1 / FL
            A = BETA * X + 4.7 * X1 / FL
            EXIT
          ENDDO
!
! STRONGLY STABLE-
!
          IF(flag_30) THEN
            IF (FL .LE. Z0) FL = Z0 + 1.E-5
            DO I = 1, 5
              ARG1 = FL / Z0
              X1 = LOG(ARG1)
# if defined (DOUBLE_PRECISION)
              X2 = DLOG(H/FL)
# else
              X2 = ALOG(H/FL)
# endif
              ARG1 = 1. - 1. / ARG1
              A = .74 * X1 + 4.7 * ARG1 + 5.44 * X2
              B = X1 + 4.7 * ARG1 + 5.7 * X2
              FL = A * S / (B*B)
              IF (FL .LE. Z0) FL = Z0 + 1.E-5
              IF (FL .GT. H) FL = H
            ENDDO
          ENDIF
!
! CALCULATE USTAR AND Z0NEW
!
        ENDIF
        TSTAR = FK * DTHETA / A
        USTAR = FK * UM / B
        Z0NEW = .00459 * USTAR * USTAR
        IF (ITER .GT. 5 .AND. ABS((USTAR - UST1)/UST1) .LT. EPS) EXIT
        UST1 = USTAR
        Z0 = Z0NEW
      ENDDO
!
! IF COME HERE, TOO MANY ITERATIONS (UGH - UGH)
!
      IF(ITER .gt. 20) THEN
        WRITE (6,70)
   70   FORMAT ('0TOO MANY ITERATIONS ON Z0 IN SUBROUTINE UZL - CHECK ','METEOROLOGICAL DATA - PROGRAM TERMINATED')
        WRITE(6,*)'U, ZM, TD, CD, CH, Z0, FL'
        WRITE(6,*)U, ZM, TD, CD, CH, Z0, FL
        STOP
      ELSE
        Z0 = Z0NEW
        CD = (USTAR/UM) ** 2
        CH = FK * FK / (A*B)
      ENDIF
      RETURN
      END SUBROUTINE UZL
# endif
END MODULE MOD_SOLAR

